6  Trabajando con datos espaciales en R

6.1 Introducción

El objetivo de este \(lab\) es familiarizarse con el uso y representación de información espacial en \(R\). Además de las funciones estadísticas que ya hemos visto, existen en \(R\) miles de paquetes, cada uno con un propósito específico. Varios de ellos permiten interactuar con datos espaciales, tanto en formato vectorial como raster. Los principales paquetes “espaciales” son:

• rgdal
• sf
• terra
• stars
• tmap

Estas librerías (paquetes) permiten leer y gestionar objetos espaciales, ya sea en formato vectorial o ficheros raster como ASCII o cualquier otro formato soportado por GDAL. Estos paquetes funcionan como un GIS, permitiéndonos realizar la mayoría de geoprocesos (intersecciones, uniones, etc.), así como procesos geoestadísticos avanzados y crear mapas.

Durante años, tres de los paquetes más populares para trabajar con datos espaciales en R han sido rgdal, raster y sp, que se convirtieron prácticamente en un estandard. Sin embargo, por diversos problemas, en 2022 se decidió deprecarlos, es decir, dejar de darles soporte, y se recomendó sustituirlos por sus “equivalentes”: terray sf. Desde octubre de 2023, aún pueden usarse las funciones de sp y raster, pero rgdal ha sido retirado, por lo que nos centraremos en conocer las funcionalidades de terray sf.

Obviamente, no podremos profundizar en toda la potencialidad de estos dos paquetes, que cada vez más logran sustituir perfectamente a un GIS de escritorio, pero sí veremos algunas de sus opciones más básicas.

6.2 Trabajando con vectores

6.2.1 Cargar información vectorial

Sabemos que los archivos shapefile de ESRI, que son casi un estándar en el mundo del SIG, son en realidad un conjunto de varios ficheros. Pero para cargarlos en \(R\), al igual que hacemos en ArcMap, basta con cargar el que tiene extensión .shp, ya que el resto van asociados a él. Sin embargo, si queremos copiar o cortar nuestros “shapes” en otra carpeta debemos tener cuidado de copiar todos los ficheros.

Veamos un ejemplo, importando las capas estaciones_meteo.shp y provincias.shp, que se encuentran en la carpeta datos/shapes, dentro de la sección “Recursos” del campus virtual de la asignatura. Estas capas contiene la información espacial de parte de las estaciones meteorológicas de las que extrajimos la información para ajustar la regresión lineal, así como los límites de todas las provincias de España.

Para leerlos en R podemos usar la función st_read() del paquete sf, por lo que antes que nada, instalamos y cargamos el paquete sf:

# install.packages('sf')
library(sf)

Y posteriormente cargamos los ficheros, especificando la ruta al archivo .shp y asignándoles un nombre:

estaciones <- st_read('data/meteo/meteo_espacial/estaciones_meteo.shp') 

Reading layer `estaciones_meteo' from data source 
  `C:\Users\Usuari\OneDrive - udl.cat\Teaching\EFCN_Geoestadística\Labs_Geoestadistica_Forestal\data\meteo\meteo_espacial\estaciones_meteo.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 90 features and 13 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 574966 ymin: 4514365 xmax: 818327 ymax: 4635648
Projected CRS: ETRS89 / UTM zone 30N

Vemos que, al leer los ficheros, la consola nos informa sobre el nuevo objeto que estamos creando. Identifica que el fichero de origen es un shapefile de ESRI, y nos informa de que se convertirá en un objeto “simple feature” con 90 observaciones de tipo “punto” y 13 campos. También nos informa sobre el bounding box del objeto espacial (el rectángulo que lo contiene) y su sistema de referencia de coordenadas o CRS.

Un objeto simple feature es el método que ha elegido sf para lidiar con información espacial. Los objetos espaciales almacenan tanto información relativa a las características espaciales (sistema de referencia, extensión, coordenadas de los objetos…) como a los atributos de cada uno de los objetos (en este caso, la información meteorológica). La ventaja de sf es que lo hace de manera muy lógica y sencilla. Vemos en el panel de Environment que estaciones aparece como un data frame normal. Y es que de hecho es un data frame normal, pero ahora es a la vez más cosas:

class(estaciones)

[1] "sf"         "data.frame"

La diferencia es que contiene, además de los atributos que veríamos en ArcMap, una columna llamada geometry con las coordenadas de cada una de las observaciones. La ventaja de esto es que el objeto espacial es un data frame, y por lo tanto podemos procesarlo con normalidad, filtrando, creando nuevas variables o modificando valores, igual que hacemos con los data frames. Al cargar estaciones nos informaba de que se traba de un objeto espacial de puntos. Por tanto cada punto contendrá las coordenadas x e y de ese punto según el CRS del objeto espacial:

estaciones

Simple feature collection with 90 features and 13 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 574966 ymin: 4514365 xmax: 818327 ymax: 4635648
Projected CRS: ETRS89 / UTM zone 30N
First 10 features:
   OBJECTID FID_weathe INDICATIVO AÃ_O MES                    NOMBRE ALTITUD
1         1         12      9382E 2012   6         PANCRUDO (D.G.A.)    1285
2         2         85      9987N 2012   6   DELTEBRE (PARC NATURAL)       4
3         3         11      9377E 2012   6               EL PEDREGAL     985
4         4         50      9567U 2012   6           EJULVE (D.G.A.)    1095
5         5         39      9530V 2012   6         UTRILLAS (D.G.A.)     960
6         6          0       3013 2012   6          MOLINA DE ARAGON    1063
7         7         70      9939A 2012   6        FUENTESPALDA (DGA)     720
8         8         40      9531X 2012   6    MONTALBAN 'AUTOMATICA'     885
9         9         84      9981A 2012   6 TORTOSA (OBSER. DEL EBRO)      48
10       10         89       9999 2012   6                      ODON    1110
   T_MAX_abs T_MIN_abs TMed_MAX TMed_MIN Tmed_MES   Provincia
1       29.5       3.0     22.0      9.6    15.80      Teruel
2       32.0      10.0     27.5     17.4    22.45   Tarragona
3       31.0       3.5     21.5     10.2    15.85 Guadalajara
4       30.0       5.0     21.9     11.5    16.70      Teruel
5       32.5       3.5     23.7     10.8    17.25      Teruel
6       32.8       3.2     23.8      8.4    16.10 Guadalajara
7       32.0       5.0     24.4     13.2    18.80      Teruel
8       31.2       5.4     23.0     11.0    17.00      Teruel
9       34.4      13.0     29.3     17.1    23.20   Tarragona
10      32.0       4.0     22.4     10.2    16.30      Teruel
                 geometry
1  POINT (666121 4514365)
2  POINT (814489 4514857)
3  POINT (620767 4515276)
4  POINT (705449 4516865)
5  POINT (681326 4520030)
6  POINT (593975 4522166)
7  POINT (758868 4522277)
8  POINT (686217 4522281)
9  POINT (794466 4524785)
10 POINT (620133 4526863)

Vamos a cargar ahora la capa provincias:

provincias <- st_read('data/meteo/meteo_espacial/provincias_spain.shp') 

Reading layer `provincias_spain' from data source 
  `C:\Users\Usuari\OneDrive - udl.cat\Teaching\EFCN_Geoestadística\Labs_Geoestadistica_Forestal\data\meteo\meteo_espacial\provincias_spain.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 51 features and 12 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -14129.47 ymin: 3892590 xmax: 1126923 ymax: 4859517
Projected CRS: ETRS89 / UTM zone 30N

El proceso es idéntico, pero ahora nos dice que el objeto creado contiene una geometría tipo multipolygon, es decir, que es un fichero vectorial de polígonos. Ahora, por lo tanto, el campo geometry contendrá, para cada observación, el conjunto de coordenadas que define los vértices del polígono:

provincias

Simple feature collection with 51 features and 12 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -14129.47 ymin: 3892590 xmax: 1126923 ymax: 4859517
Projected CRS: ETRS89 / UTM zone 30N
First 10 features:
   OBJECTID                INSPIREID COUNTRY
1         1 ES.IGN.BDDAE.34160100000      ES
2         2 ES.IGN.BDDAE.34080200000      ES
3         3 ES.IGN.BDDAE.34100300000      ES
4         4 ES.IGN.BDDAE.34010400000      ES
5         5 ES.IGN.BDDAE.34070500000      ES
6         6 ES.IGN.BDDAE.34110600000      ES
7         7 ES.IGN.BDDAE.34040700000      ES
8         8 ES.IGN.BDDAE.34090800000      ES
9         9 ES.IGN.BDDAE.34070900000      ES
10       10 ES.IGN.BDDAE.34111000000      ES
                                                                        NATLEV
1  https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
2  https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
3  https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
4  https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
5  https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
6  https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
7  https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
8  https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
9  https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
10 https://inspire.ec.europa.eu/codelist/AdministrativeHierarchyLevel/3rdOrder
   NATLEVNAME     NATCODE         NAMEUNIT CODNUT1 CODNUT2 CODNUT3 Shape_Leng
1   Provincia 34160100000      Araba/Álava     ES2    ES21    <NA>   639475.2
2   Provincia 34080200000         Albacete     ES4    ES42    <NA>   770993.4
3   Provincia 34100300000 Alacant/Alicante     ES5    ES52    <NA>   583497.3
4   Provincia 34010400000          Almería     ES6    ES61    <NA>   632690.0
5   Provincia 34070500000            Ávila     ES4    ES41    <NA>   630836.7
6   Provincia 34110600000          Badajoz     ES4    ES43    <NA>  1136844.0
7   Provincia 34040700000    Illes Balears     ES5    ES53    <NA>  1459098.3
8   Provincia 34090800000        Barcelona     ES5    ES51    <NA>   850790.3
9   Provincia 34070900000           Burgos     ES4    ES41    <NA>  1138046.5
10  Provincia 34111000000          Cáceres     ES4    ES43    <NA>   965487.3
    Shape_Area                       geometry
1   3035458537 MULTIPOLYGON (((519021.3 47...
2  14920858305 MULTIPOLYGON (((539277 4215...
3   5820590272 MULTIPOLYGON (((697663 4195...
4   8767315137 MULTIPOLYGON (((496730.8 39...
5   8047820788 MULTIPOLYGON (((292983.8 44...
6  21793269761 MULTIPOLYGON (((165751.6 42...
7   5018967692 MULTIPOLYGON (((868245.7 43...
8   7762160343 MULTIPOLYGON (((898327 4573...
9  14279022189 MULTIPOLYGON (((419298.4 46...
10 19885694014 MULTIPOLYGON (((172175.6 43...

provincias$geometry[1]

Geometry set for 1 feature 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 476582.6 ymin: 4702304 xmax: 562700.7 ymax: 4784934
Projected CRS: ETRS89 / UTM zone 30N

MULTIPOLYGON (((519021.3 4717987, 518976.3 4717...

6.2.2 Representando gráficamente información vectorial

Evidentemente, cuando trabajamos con datos espaciales una de las acciones más interesantes es representarlos gráficamente. Para ello basta con usar la función plot() y el nombre del objeto espacial

plot(estaciones)

Warning: plotting the first 9 out of 13 attributes; use max.plot = 13 to plot
all

Por defecto nos ploteará todos los campos. Si queremos visualizar específicamente alguno de ellos debemos indicarlo expresamente, poniendo el nombre del campo entre corchetes []:

plot(estaciones["Tmed_MES"])

\(R\) mapea los puntos según su ubicación, y les asigna un color en función de los valores de la variable elegida, pero podemos editar la visualización igual que hacemos con los \(scatterplot\) o los \(plots\) normales.

Por ejemplo, podemos cambiar el tipo de símbolo con pch, el color con col o el tamaño con cex.

plot(estaciones["Tmed_MES"], pch = 19)

plot(estaciones["Tmed_MES"], pch = 21, col='black', bg= 'red')

En caso de que sólo queremos plotear la geometría del objeto (sin incluir ningún atributo) debemos indicarlo con la función st_geometry()

plot(st_geometry(provincias))

plot(st_geometry(provincias), col = "dark red")

Podemos también combinar dos capas diferentes mediante el comando add = TRUE:

plot(st_geometry(provincias) )
plot(estaciones["Tmed_MES"], pch = 19, col = "red", cex = 0.5, add = TRUE) 

6.2.2.1 Visualizando capas vectoriales con tmap

Un paquete interesante para la visualización de información espacial es tmap. Este paquete funciona de manera diferente, ya que le debemos indicar las diferentes capas a visualizar uniendo las órdenes con el símbolo +. Lo primero es indicar la capa que queremos visualizar mediante el comando tm_shape() y luego especificamos si queremos ver puntos (tm_dots()), los bordes de los polígonos (tm_borders()), o los polígonos con relleno (tm_polygons(), entre otras muchas opciones). De cada uno de los comandos tm_* puede personalizarse el color, forma, tamaño, etc. con valores fijos o en función de alguna columna.

Probemos a visualizar las estaciones:

library(tmap)

Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
remotes::install_github('r-tmap/tmap')

tm_shape(estaciones) +
    tm_dots()

Lógicamente, podemos cambiar el tamaño de los puntos, su forma, o los colores, haciéndolos depender de la variable Tmed_MES e incluso definir la paleta:

tm_shape(estaciones) +
    tm_dots(size = 1, shape = 20, col = "Tmed_MES", palette = "viridis")

Con tmaps podemos visualizar tantas capas como queramos, sólo necesitamos volver a utilizar la función tm_shape() con el nombre de la nueva capa. Eso sí, la extensión del mapa vendrá dada por la primera capa que llamemos. Probemos a visualizar las estaciones sobre las provincias del nordeste:

tm_shape(provincias) +
    tm_polygons(col = "lightgrey") +
    tm_text("NAMEUNIT") +
tm_shape(estaciones) +
    tm_dots(col = "Tmed_MES", palette = "viridis", size = 0.75)

La extensión espacial del mapa vendrá determinada por la primera capa o shape que carguemos. Si cargamos primero las estaciones:

tm_shape(estaciones) +
    tm_dots(col = "Tmed_MES", palette = "viridis", size = 0.75) +
tm_shape(provincias) +
    tm_borders() +
    tm_text("NAMEUNIT") 

Las opciones generales del mapa se pueden personalizar mediante tm_layout():

tm_shape(estaciones) +
    tm_dots(col = "Tmed_MES", palette = "viridis", size = 0.75,title = "T media") +
tm_shape(provincias) +
    tm_borders() +
    tm_text("NAMEUNIT") +
    tm_layout(legend.outside = F, bg.color = "steelblue",title = "Estaciones meteorológicas del Nordeste", title.size = 4)

Incluso podemos definir un fondo basado en un proveedor como Google Maps o OpenStreetMap, con la función tm_basemap() y hacer que el mapa sea interactivo definiendo antes tmap_mode("view"): (una lista de las opciones se puede encontrar aquí)

tmap_mode("view")

tmap mode set to interactive viewing

tm_shape(estaciones) +
    tm_dots(col = "Tmed_MES", palette = "viridis") +
tm_shape(provincias) +
    tm_borders() +
    tm_text("NAMEUNIT") +
    tm_basemap("OpenStreetMap.Mapnik")

Por último, una vez tengamos un mapa que nos guste, podemos exportarlo al formato de nuestra elección:

mapa <- tm_shape(estaciones) +
    tm_dots(col = "Tmed_MES", palette = "viridis") +
    tm_shape(provincias) +
    tm_borders() +
    tm_text("NAMEUNIT") +
    tm_basemap("OpenStreetMap.Mapnik")

## Guardar como imagen (modo "plot")
tmap_save(mapa, filename = "mapa_estaciones.png")

## Guardar como HTML file (modo "view")
tmap_save(mapa, filename = "mapa_estaciones.html")

6.2.2.1.1 Para saber más

Existen numerosos tutoriales online con las principales funcionalidades del paquete tmap. Algunas de las recomendadas son tmap: get started!, el libro tmap book o el capítulo dedicado a tmap del libro Making maps with R

6.2.3 Seleccionando campos, filtros y observaciones

Pero ¿qué pasa si no queremos trabajar con todos los puntos, o si solo queremos representar una parte de los mismos? Como comentábamos antes, un objeto sf no es más que un data.frame con geometría espacial, por lo que podemos hacer con el todas las operaciones que hacemos con un data frame. Por ejemplo, podemos escoger submuestras basadas en valores determinados mediante la función subset():

bajas <- subset(estaciones, ALTITUD <= 500)

tm_shape(estaciones) +
    tm_dots(col = "darkgrey") +
tm_shape(bajas) +
    tm_dots(col= "Tmed_MES")

Este comando lo podemos usar también con caracteres de texto:

catalunya <- subset(provincias, NAMEUNIT %in% c("Lleida", "Tarragona"))

tm_shape(catalunya) +
    tm_borders(col = "orange") +
tm_shape(estaciones) +
    tm_dots(col= "Tmed_MES")

Sin embargo, hay que tener en cuenta que si lo que hacemos es seleccionar una columna extraeremos los valores de la variable como un vector, pero perderemos la información espacial:

estaciones$Tmed_MES

 [1] 15.80 22.45 15.85 16.70 17.25 16.10 18.80 17.00 23.20 16.30 17.45 18.75
[13] 18.45 16.25 21.75 19.75 17.30 18.75 16.10 19.75 19.00 17.90 19.90 20.35
[25] 16.50 19.90 20.70 22.75 16.35 18.85 16.15 21.65 20.65 20.80 20.80 22.30
[37] 19.60 22.05 16.40 22.95 20.60 18.70 18.30 19.60 23.50 23.35 22.15 17.80
[49] 18.35 21.30 16.80 18.00 17.40 18.00 21.60 19.15 19.35 19.55 21.75 20.85
[61] 21.00 19.85 20.15 20.90 21.70 18.80 16.95 18.10 22.15 22.15 20.55 22.10
[73] 21.90 21.75 21.95 20.90 21.00 20.55 21.20 20.70 20.30 19.80 20.00 19.85
[85] 20.70 20.80 20.70 16.75 22.05 21.50

6.2.4 Operaciones espaciales con sf

sf también permite realizar todas las operaciones espaciales típicas de un GIS como intersectar capas, unirlas, crear buffers, etc. No vamos a ver todo esto con mucho detalle porque se escapa del objetivo de la asignatura, pero podéis encontrar un tutorial muy extenso aquí.

De momento solo veremos un ejemplo de intersección entre dos capas:

estaciones_cat <- st_intersection(catalunya, estaciones)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

tm_shape(catalunya) +
    tm_borders(col = "orange") +
tm_shape(estaciones_cat) +
    tm_dots(col= "Tmed_MES")

6.3 Trabajando con rasters

La información en formato raster es básicamente una matriz de dos dimensiones, en la que cada una de las celdas o píxeles tiene un valor numérico que representa una variable, numérica o categórica. Este tipo de estructura es muy adecuada para representar información sobre fenómenos continuos, como la temperatura, elevación, distancias…

6.3.1 Cargando rasters

Igual que hicimos con las capas vectoriales, lo primero de todo es cargar la capa en un objeto de R. En este caso usaremos la funcion rast(), del paquete terra. Vamos a cargar un modelo digital de elevaciones de la zona de estudio de las estaciones meteorológicas que se encuentra en la carpeta “Recursos/data/rasters_meteo” del campus virtual

# install.packages('terra', dep = TRUE)

library(terra)

terra 1.7.78

dem <- rast("data/meteo/meteo_espacial/rasters_meteo/elevation.txt")
dem

class       : SpatRaster 
dimensions  : 3598, 2594, 1  (nrow, ncol, nlyr)
resolution  : 100, 100  (x, y)
extent      : 561200, 820600, 4398400, 4758200  (xmin, xmax, ymin, ymax)
coord. ref. :  
source      : elevation.txt 
name        : elevation 

Y podemos visualizarlo de forma sencilla mediante plot():

plot(dem)

6.3.2 Cargar varias capas raster

Puede que os hayáis dado cuenta de que la carpeta rasters_meteo contiene más ficheros además del de elevaciones. Gracias a las utilidades del paquete terra podremos cargar varias capas raster a la vez, que quearán almacenadas en lo que se llama stack que no es otra cosa que un objeto que contiene numerosas capas raster que comparten extensión, resolución, sistema de coordenadas… Un ejemplo típico serían las distintas bandas de una imagen de satélite.

Sólo necesitamos una lista que contenga los nombres de los objetos a cargar (rasters). Podemos usar para ello la función list.files():

lista <- list.files('data/meteo/meteo_espacial/rasters_meteo/',full.names = TRUE)
lista

[1] "data/meteo/meteo_espacial/rasters_meteo/d_atl.txt"    
[2] "data/meteo/meteo_espacial/rasters_meteo/d_medit.txt"  
[3] "data/meteo/meteo_espacial/rasters_meteo/elevation.txt"
[4] "data/meteo/meteo_espacial/rasters_meteo/lat.txt"      
[5] "data/meteo/meteo_espacial/rasters_meteo/long.txt"     

list.files devuelve un vector con los nombres de todos los ficheros dentro de una determinda carpeta. Esta función tiene algunos argumentos interesantes:

• full.names: TRUE o FALSE determina si queremos que devuelva sólo el nombre del archivo o mejor la ruta completa a cada uno de ellos.
• pattern: parámetro que permite filtrar los objetos por nombre. Por ejemplo, si queremos sólo los ficheros en formato txt:

lista <- list.files('data/meteo/meteo_espacial/rasters_meteo/',full.names = TRUE, pattern = 'txt')
lista

[1] "data/meteo/meteo_espacial/rasters_meteo/d_atl.txt"    
[2] "data/meteo/meteo_espacial/rasters_meteo/d_medit.txt"  
[3] "data/meteo/meteo_espacial/rasters_meteo/elevation.txt"
[4] "data/meteo/meteo_espacial/rasters_meteo/lat.txt"      
[5] "data/meteo/meteo_espacial/rasters_meteo/long.txt"     

Una vez que tenemos una lista de los rasters a cambiar, podemos cargar la lista con la función rast():

rasters <- rast(lista)

Y si ploteammos veremos que hay varios raster (capas) cargadas:

plot(rasters)

Podemos acceder a cada una de las “bandas” del stack con el comando $:

plot(rasters$elevation)

6.4 Proyección y sistema de referencia

El sistema de referencia (CRS) es un elemento clave de la información espacial. Hasta ahora hemos trabajado la información espacial sin prestar atención a este parámetro. Sin embargo, tarde o temprano tendremos que ocuparnos de él. No tenerlo en cuenta nos dará problemas. Por ejemplo, no podremos superponer o combinar capas por la ausencia o diferencia de sistemas de referencia, por lo que tenemos que saber cómo ensamblar o reproyectar una capa a un sistema de referencia diferente.

6.4.1 Asignar un sistema de referencia de coordenadas (CRS)

Un CRS sólo deber asignarse a una capa cuando esta carezca de esta información de manera explícita. Lo primero es comprobar por tanto qué CRS tienen las capas con las que estamos trabajando:

crs(rasters)

[1] ""

crs(estaciones)

[1] "PROJCRS[\"ETRS89 / UTM zone 30N\",\n    BASEGEOGCRS[\"ETRS89\",\n        ENSEMBLE[\"European Terrestrial Reference System 1989 ensemble\",\n            MEMBER[\"European Terrestrial Reference Frame 1989\"],\n            MEMBER[\"European Terrestrial Reference Frame 1990\"],\n            MEMBER[\"European Terrestrial Reference Frame 1991\"],\n            MEMBER[\"European Terrestrial Reference Frame 1992\"],\n            MEMBER[\"European Terrestrial Reference Frame 1993\"],\n            MEMBER[\"European Terrestrial Reference Frame 1994\"],\n            MEMBER[\"European Terrestrial Reference Frame 1996\"],\n            MEMBER[\"European Terrestrial Reference Frame 1997\"],\n            MEMBER[\"European Terrestrial Reference Frame 2000\"],\n            MEMBER[\"European Terrestrial Reference Frame 2005\"],\n            MEMBER[\"European Terrestrial Reference Frame 2014\"],\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]],\n            ENSEMBLEACCURACY[0.1]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4258]],\n    CONVERSION[\"UTM zone 30N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-3,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"Europe between 6°W and 0°W: Faroe Islands offshore; Ireland - offshore; Jan Mayen - offshore; Norway including Svalbard - offshore; Spain - onshore and offshore.\"],\n        BBOX[35.26,-6,80.49,0.01]],\n    ID[\"EPSG\",25830]]"

En este caso, los rasters que hemos cargado no tienen asignado ningún CRS, así que debemos asignarle uno. Para ello podemos comprobar la lista de los diferentes CRS - y manera de codificarlo - en esta web. EN este caso sabemos que el CRS correcto es EPSG: 23030 - UTM ED50 30N, que en el codificado “proj4” es:

+proj=utm +zone=30 +ellps=intl +units=m +no_defs

Y lo asignamos como:

crs(rasters)<- "+proj=utm +zone=30 +ellps=intl +units=m +no_defs"

Ahora podemos comprobar CRS otra vez:

crs(rasters)

[1] "PROJCRS[\"unknown\",\n    BASEGEOGCRS[\"unknown\",\n        DATUM[\"Unknown based on International 1924 (Hayford 1909, 1910) ellipsoid\",\n            ELLIPSOID[\"International 1924 (Hayford 1909, 1910)\",6378388,297,\n                LENGTHUNIT[\"metre\",1,\n                    ID[\"EPSG\",9001]]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8901]]],\n    CONVERSION[\"UTM zone 30N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-3,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]],\n        ID[\"EPSG\",16030]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]]]"

6.4.2 Proyectar a otro CRS

Es común que, cuando descargamos cartografía de diversas fuentes, cada una tenga su propio CRS. Esto nos trae problemas en cuanto a la visualización, pero también si queremos hacer operaciones espaciales (intersects, buffers…). Para solventarlo, debemos proyectar alguna de las capas, de manera que estén todas en la misma CRS.

NOTA: es importante tener en cuenta que proyectar una capa a un nuevo CRS no es lo mismo que asignar un CRS. En el primer caso, la capa ya tiene CRS asignado, y lo que hacemos es transformar espacialmente las coordenadas. En el segundo (asignar) las coordenadas están, pero no está definido que CRS usar para representar.

nuevo_provincias <- st_transform(provincias, "+proj=utm +zone=31 +ellps=intl +units=m +no_defs")

plot(st_geometry(nuevo_provincias))
plot(estaciones["Tmed_MES"], pch = 19, col = "red", cex = 0.5, add = TRUE) 

Sin embargo, el mismo mapa con tmaps sí que funcionará:

tmap_mode("plot")

tmap mode set to plotting

tm_shape(nuevo_provincias)+
    tm_borders()+
    tm_shape(estaciones)+
    tm_dots()

Esto es porque tmaps implementa, igual que ArcMap la llamada proyección on-the-fly, es decir que transforma automáticamente los CRS de las capas al CRS definido al principio (fijaos que las provincias salen ahora giradas respecto a la visualización original).

6.5 Visualización de mapas combinados de rasters y shapefiles

Al tener asignado correctamente el CRS; ahora podemos plotear las diferentes capas a la vez usando las funciones de tmap:

# tmap_mode("plot")

tm_shape(rasters$elevation) +
    tm_raster() +
tm_shape(provincias) +
    tm_borders() +
tm_shape(estaciones) + 
    tm_dots( col="Tmed_MES", size = 0.5, palette= "viridis") +
    tm_layout(legend.outside = T)

6.6 Para saber más

Desde luego, esto sólo es una introducción al uso de \(R\) como entorno de procesado y visualización SIG. Hay mucho más, y para profundizar os recomiendo consultar los siguientes recursos:

• Visualización con tmap:
  • Elegant and informative maps with tmap
  • tmap: get started!

web de R Spatial (https://rspatial.org/#google_vignette) o la web Introduction to GIS with R